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Abstract 

The curse of dimensionality (COD) limits the current state-of-the-art ab initio propagation 
methods for non-relativistic quantum mechanics to relatively few particles. For stationary structure 
calculations, the coupled-cluster (CC) method overcomes the COD in the sense that the method 
scales polynomially with the number of particles while still being size-consistent and extensive. We 
generalize the CC method to the time domain while allowing the single-particle functions to vary in 
an adaptive fashion as well, thereby creating a highly flexible, polynomially scaling approximation 
to the time-dependent Schrbdinger equation. The method inherits size-consistency and extensivity 
from the CC method. The method is dubbed orbital-adaptive time-dependent coupled-cluster 
(OATDCC), and is a hierarchy of approximations to the now standard multi-configurational time- 
dependent Hartree method for fermions. A numerical experiment is also given. 
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I. INTRODUCTION 



Presently, the most advanced ab initio approximations to the time-dependent Schrodinger 
equation for a system of identical particles are the multiconfigurational time-dependent 
Hartree methods for fermions (MCTDHF) and variants l|-|3|]. These methods apply the 
time-dependent variational principle {4-6] to an iV-body wavefunction ansatz being a Slater 
determinant expansion using a finite (incomplete) set of L > N orbitals tp p with creation 
operators cj, (with {c p ,cj} = 5 pq for fermions), 

I^mctdhf) = E E " ' E A pi-pn c IA2 ■ ■ - C LH' 

Pi P2>Pl PN>PN-1 

where both the amplitudes A pi ... PN and the orbitals <p p are free to vary in single-particle 
space. Varying the orbitals in this way is especially important if studies of unbound systems 
are desired, such as the study of ionization of atoms or molecules. The key point is that, 
if the orbitals are not optimized, a system in the continuum would need a huge fixed basis. 
Moreover, the variational determination of the orbitals compresses the wavefunction in a 
quasi-optimal way: from time t to t + alt, the basis is changed as so to minimize the L 2 norm 
error of the wavefunction 6]. 

While powerful, MCTDHF still suffers from exponential scaling of computational com- 
plexity with respect to the number of particles N present; the somewhat prosaically termed 
"curse of dimensionality" (COD). The effect of the variational determination of the single- 
particle functions can be said to be a postponing of the COD to higher particle numbers. 
For example, for the simple time-dependent Hartree-Fock method (TDHF), i.e., MCTDHF 
using precisely L = N orbitals and therefore only a single determinant \<f>), qualitatively 
good results may be achieved, even if the system is unbound 3). 

One may attempt at reducing the exponential scaling by truncating the Slater determi- 
nant expansion at, say, single and double excitations relative to one of the determinants 
\4>) = c\ ■ • • c* N \ — ), considered as a "reference determinant", 



I* 



MCTDHF-SD 



) = (i + E^ c ^) + iE< c ^ c M |0) = (i+a sd )|0), (i) 



2! 

ijab 



hoping that the higher-order excited determinants' contribution can be neglected. (We have 
arbitrarily chosen intermediate normalization (<f)\^) = 1. In the sums, i,j < N and a,b > 
N is assumed.) This would achieve polynomial scaling but would destroy the important 



property of size-consistency js, 9|: approximation of non-interacting subsystems separately 
at the singles and doubles level would not be consistent with approximating the whole at the 
same level; independent excitations of the subsystems are neglected, giving rise to artificial 
correlation effects. 

In this article, we develop a time- dependent version of the popular coupled-cluster (CC) 
method for fermions, where we allow the orbitals to vary in a similar fashion to MCTDHF. 
We call the method orbital adaptive time- dependent coupled- cluster (OATDCC). Formally, 
the two methods are very similar, with closely related equations of motion. The main 
difference is the fact that coupled-cluster is not variational in the usual sense, rather, it is 



naturally cast in a bivariational setting, a generalization of the variational approach [10]. 
Bivariational functionals are complex analytic, while the standard variational functionals 
are manifestly real. Moreover, approximations of both the wavefunction |\&) the complex 
conjugate must be introduced. For the version of CC that has now become standard 
(referred to as "standard CC" in this paper), the wavefunctions in the singles and double 
approximation (CCSD) are parametrized according to 



|* cc >^e T |0>, T = Y j rtc^ l + ^Y,r^c\c4c 3 + --- (2a) 

ai ijab 

(feel = (01(1 + A)e" T , A = A a*a + ^ b c\c a c]c b + ■■■ (2b) 

ai ijab 

where (\&cc| approximates (\]/|/(\]/|\l/). The operator T is called a cluster operator and 
produces excitations with respect to a reference Slater determinant The amplitudes rf 
and Tj" 6 correspond (to first order) to the expansion coefficients A® and of Eqn. (|T|). The 
cluster operator A is a de-excitation operator, and its amplitudes and A^ 6 are essentially 
the parameters of (^ccIj which also is composed of excitations relative to a reference bra 
determinant (<j>\. To anticipate the developments in later sections, we have introduced 
creation and annihilation operators with respect to biorthogonal orbitals <p p and (p q , i.e., 

{c p , c\) = c p c\ + c\c p = (<p p \(p q ), 

That is to say, |\&) is built using the (f p , while is built using (p v . This relaxation of 
orthonormality of the orbitals is necessary to ensure that the bivariational functional is 
complex analytic if the orbitals are to be treated as variational parameters, as discussed in 
Section [mi 
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Our treatment can be viewed as a generalization of the standard CC Lagrangian approach 



to linear response theory 



10Ml3j. where the A amplitudes are time- dependent Lagrangian 



multipliers as introduced in a constrained minimization of the CC energy. However, we 
emphasize the bivariational point of view, where A becomes a part of the wave function 
parametrization. In Ref. [l^j even biorthogonal orbitals were consideres. 

For an excellent introduction to CC theory, see the article 15] by Crawford, and the 
review 16] by Bartlett and Musial, as well as the textbooks {9] by Shavitt and Bartlett and 



[17] by Harris et al. For the present work, the article [10] is fundamental, as it casts the CC 
theory in the bivariational framework, an approach not emphasized in most introductions 
to CC theory. 

The OATDCC method considers the standard CC ansatz plugged into the proper bivari- 
ational functional, see Section fill Dl below. In addition to having the T and A amplitudes as 
degrees of freedom, the orbitals \ip p ) and approximations (<p p \ to their complex conjugates 
are varied. It turns out the relaxation of the orbitals makes the singles amplitudes rf and \\ 
redundant. Truncation of the remaining terms at doubles, triples, etc, then gives a hierarchy 
of approximations with TDHF at one end (T = A = 0), and full MCTDHF (no truncation of 
T or A) at the other. Inbetween we have the doubles approximation (OATDCCD), doubles- 
and-triples approximation (OATDCCDT), and so on. Considering the success of the CC 
method for structure calculations and the MCTDHF method for dynamics, OATDCC should 
be a viable alternative to MCTDHF with asymptotically much lower cost but good accuracy. 
Importantly, OATDCC is size-consistent. 

There are few applications of CC methods to ab initio dynamics in the ^literature. It 
was, however, proposed as early as 1978 by Schonhammer and Gunnarsson 



and inde- 



pendently by Hoodbhoy and Negele 



19 



20] . who even considered time-dependent orbitals, 



albeit with an explicit time-dependence. We shall discuss their approach briefly in Section 
IV Al Recently, the standard CC ansatz using a fixed basis was applied to laser-driven dy- 
namics of some small molecules 21], but with expectation values calculated in a different 
way from the usual CC approach. 

This article is written with the MCTDHF community in mind. Since bivariational princi- 
ples are rarely considered (which is true for the CC community as well), we give a somewhat 
detailed discussion in Section [Til As CC theory might be unfamiliar, and as we use a some- 
what unfamiliar variational approach to derive CC theory, Section [TTTJ is devoted to the 



basics of the CC formalism, leading up to the case where the orbitals are varied freely. In 
Section IIVI we discuss the OATDCC functional and derive the corresponding equations of 
motion. We then perform a simple numerical experiment in Section I VI I to demonstrate the 
method before we conclude the paper. 

CC calculations invariably involve a lot of algebra. Therefore, in |XJ algebraic expressions 
for various quantities appearing in the OATDCCD method are listed. These are generated 
using symbolic algebra software developed with the SymPy library for the programming 



language Python [22]. An independent derivation of MCTDHF is given in|B]in order to 
shed further light on the connections between the variational and the bivariational principles. 
It may also serve as a helpful device for the researchers in the audience not familiar with 
MCTDHF theory. 

No attem pt is m ade to be mathematically rigorous in this article, as like the standard 



CC method 
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251 ] and the MCTDHF method [y, |26H28| such an analysis is expected to 
be quite involved. Instead, we make formal computations as if all operators present were 
bounded or the spaces finite dimensional. 

II. BIVARIATIONAL PRINCIPLES 
A. Functionals 

Let A be an operator over Hilbert space H, and consider the functional 

defined whenever the expression makes sense. Note that the arguments are two independent 
Hilbert space elements. The functional E-a is a generalization of the expectation value 
functional to operators that are not necessarily Hermitian, and for obvious reasons it may 
be called the bivariational expectation value functional fiol . I29I ]. Consider the conditions 
for vanishing first variation, SSa = 0, for all independent variations of (\&'| and \^/). A 
straightforward formal calculation gives the stationary conditions 

(A-a)\V) = and (*'|(A-a) = 0, (3) 

with 

a = £i[<m*>] 



being the value of E-a at the critical point. In other words and |\&) with 7^ are 

left- and right eigenvectors, respectively, of the operator A, with eigenvalue a. Computing 
the eigenvalues and eigenvectors from 88a = is therefore referred to as "the bivariational 
principle." 

Since the system Hamiltonian H = W is the generator for the time-evolution of the 
system, we consider the following bivariational generalization of the usual action functional 
i.e., 

r , „ , „ f T W(t)\ (ihg-H) \m)) 

= jf ul S©- £ - tWt)l ' l * (t » ldt ' (4) 

where it is understood that the functional depends on the whole history of the system from 
time t — to t — T. Suppose S is stationary (5S = 0) under all variations of (^'| and 
|\&) vanishing at the endpoints t — 0, T. Straightforward manipulations now give, up to 
irrelevant time-dependent phase constants, 

ih^(t)) = H\V(t)) and -ih^'(t)\ = (^(t)\H. (5) 

Consequently, we note that the time- dependent Schrodinger equation and its complex con- 
jugate arise from a time- dependent bivariational principle. 

In both the stationary and time-dependent case, it is convenient to do a reparametrization 
(§| = (^'l^)" 1 ^'! so that 

<*|*> = 1, (6) 
eliminating the denominator in each functional. For the stationary case, this effectively 
eliminates one of the two eigenvalue equations ([3]), and makes a unique function of 
|^) and vice versa. To see this, note that since (a) eigenvectors corresponding to different 
eigenvalues are always orthogonal, and since (b) (^|^) = 1, (^| is uniquely given by at 
the critical point as the biorthogonal left eigenvector corresponding to the right eigenvector 

For the time- dependent case, the normalization (&(t)\ty(t)) = 1 eliminates one of the 
Schrodinger equations (jSJ), but the phase ambiguity is still present for the remaining equation, 
a similar situation as the stationary case. 

Unless otherwise stated, the normalization ([6]) is assumed in the following, and it is 
indicated by the tilde (W| instead of the prime (\&'|. 
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It is also natural to assume some secondary normalization on such that the critical 
point actually becomes locally unique for both the time-dependent and time-independent 
cases. ("Locally unique" means that there may be many critical points but that they are 
isolated.) This is done in CC theory, where (<p\^) = 1 is assumed. This removes phase 
ambiguity in both time-dependent and time-independent pictures. 



B. Generating approximations from submanifolds 

Contrary to the usual variational principles, in the bivariational principles the wavefunc- 
tion and its complex conjugate are formally independent. When applied to Hermitian oper- 
ators this opens up possibilities for more general approximations compared to the standard 
variational principles to spectra and dynamics, when the variations are restricted to some 
predefined approximation manifold. However, where the usual "Hermitian" time-dependent 
variational principle is, to paraphrase Kramer and Saraceno Q , a deaf and dumb procedure 
that always gives an answer, the bivariational principle requires a more careful approach, as 
we will discuss in this section. 

Perhaps for this reason, the bivariational principles are little known. The author has 



found only a 
cluster paper 



ew relevant sources in the literature apart from Arponen's seminal coupled- 



10] , the most relevant ones being a brief mention by Killingbeck in his review 



on perturbation theory [30] and a discussion by Lowdin et al. [29| concerning self-consistent 
field-calculations on non-Hermitian (complex scaled) Hamiltonians. The bivariational prin- 
ciple seems largely unexplored. 

It is important to note that the standard critique of coupled-cluster is that it is "non- 
variational" . While it is true that it makes estimation of errors harder, it is not a serious 
drawback in any other sense, since the calculation is firmly rooted in a variational principle. 
For example, if H does not depend explicitly on time, (d/dt)Sn = 0, i.e., energy is conserved. 
Probability is always conserved, {d/dt)£\ = (d/dt) (^\^) = 0. Just like the usual time- 
dependent variational principle, these are simple consequences of the symmetries of the 
action functional. 

For bivariational approximations, one introduces different parametrizations of the wave- 
function and its complex conjugate ($\, which is contrary to the usual variational 
principle. Notice that in Eqn. (T5]), different parametrizations |^cc) an d (^ccl are used, 



but (^ccl^cc) = 1- Formally, we do a variation over a manifold Ai C H' x H, i.e., 
((^|, |^)) G Ai. Alternatively, one may think of Ai as a subset of rank-one density opera- 
tors u = Tr(u) = 1, where ^ u is allowed. 

As Ji is a complex space, the functionals Sh and 5 are complex. On the other hand, in the 
usual variational principle, the functionals are always real-valued. This has some interesting 
consequences. We now briefly discuss four important aspects: the analytic structure of the 
functionals, the need for systematic refmability of Ai, interpretations of complex critical 
points, and even- dimensionality of Ai. 

The parametrization of u = ((^f\,\^f)) must be complex analytic, at least locally. A 
parametrization which is not analytic is, in essence, a real parametrization, since it depends 
on both the real and imaginary parts of the (local) coordinates z G C n separately, and not 
only Re z + i Im z. Thus, we have 2n real coordinates. Since both Re SS and Im 5S must van- 
ish, this leads to An equations. Unless there is some extra structure, i.e., that lm.S = such 
as in the standard variational principle, a solution cannot be expected to exist. Correspond- 
ingly, the bivariational functionals should be complex analytic. Nowhere should explicitly 
real parameters occur, and nowhere should parameters be explicitly complex conjugated. 

Suppose the system Hamiltonian H is bounded from below, i.e., the expectation value 
is bounded from below. This is the source of the usefulness of the Hermitian stationary 
variational principle, since any parametrization gives an upper bound for the ground state 
energy. A potential danger with the bivariational expectation value functional is that it is 
not bounded from below, even if H is. Indeed, Sh is complex analytic and can take on 
values in the whole of C. One cannot insert "just anything" and hope to get sensible results 
by computing critical points. To avoid this problem, and to allow for the computation of 
error estimates, Ai should be chosen in a way that is in some sense systematically refmable 
towards the full space %' x T-i, e.g., there is some discretization parameter that can be used 
to measure the accuracy. In CC theory, this parameter is the truncation level of the cluster 
operators and the number L of orbitals used. 

Critical values of Sh for H = W may be complex, even though the exact eigenvalues 
are always real. However, if H = H\ the imaginary values of the critical values generally 
"should be small" if Ai is chosen "well enough" , and may correspondingly be ignored in order 
to assign a physical interpretation to the critical value, i.e., energy. This is also justified by 
the fact that the functional Re Sh has the same critical points as Sh if the parametrization 
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is analytic. 

For the approximate manifold Ai, we must be certain that the critical point ((^|, |^)) 
is locally unique. In particular, this is necessary for the corresponding critical value, i.e., 
£a[(\I/||\I/)], to be unique for any observable A. Otherwise, the physical state cannot be said 
to be well-defined. Intuitively, the parameters must then come in pairs; roughly stated every 
parameter in should have a parameter in 

This can be shown explicitly. Suppose (W| and are parametrized locally using some set 
of complex variables z(t) = C n , i.e., we have an approximating manifold Ai C ~H! x % whose 
dimension is assumed to be finite for simplicity. When inserted into the time-dependent 
bivariational functional, we obtain a new functional = S[(^(z(-))\, |\&(z(-))}], whose 

stationary point is given by the solution of the differential equation 

ihC(z)z = V z E(z), E(z) = £h[(V(z)\, \V(z))]. (7) 

The matrix C(z) is given by 





dzj 



which is complex anti-symmetric. For any anti-symmetric matrix, if A is an eigenvalue of 
multiplicity m, so is —A, implying that C(z) is not invertible if n is odd, since it must 
have a zero eigenvalue. Consequently, the approximation manifold must be complex even 
dimensional if Equation ([7]) is to have a unique solution. 



III. COUPLED-CLUSTER FUNCTION ALS 



A. Biorthogonality of orbitals 

The fact that we require our bivariational functionals to be analytic is important, as it 
necessitates the relaxation of the ort honor mality the orbitals <p p by introducing a second set 
of biorthogonal orbitals (p q , being in effect approximate complex conjugates of each other. 
However, they are introduced as independent complex parameters. 

Consider a subspace V C "H, generated by a finite set of orbitals $ = (<pi,(p2, • • • ,<Pl), 
which we write V = V[$]. (We interpret (f p as the pth column of the matrix $.) These 
orbitals need not be orthonormal; it is only the one-body space spanned by the <p p that 
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matters. In the bivariational functionals, we desire to vary in as large space V C w 
as possible, while guaranteeing the existence of a dual vector non-orthogonal to |\&) £ V. 
(Otherwise, the denominator in Eqn. (jlj) may vanish.) The only restriction on the space V 
is that it is generated by a set of dual orbitals $ = (<pi, ■ ■ ■ ; <£l) (where we interpret (p v as 
the pth row of $). Sometimes we will stress the fact that ip p and (p p are single-particle ket 
and bra- functions, respectively, by explicitly writing \ip p ) and {<p p \. 

Consider the overlap matrix S with matrix elements S pq = (tp p \tf q ). Since the spaces 
V[$] and V[<3>] only depend on the subspaces spanned by $ and $, respectively, we may via 
a suitable transform assume that S is diagonal with only Is and Os on the diagonal. If S is 
invertible, then the orbitals are biorthogonal, 

(<Pp\<Pq) = <W 

It is straightforward to show the following claim: the existence of a (\&| £ V[$] for every 
|\&) £ V[$] such that (^1^) 7^ is equivalent to requiring the overlap matrix S pq = (tp p \(p q ) 
to be invertible, i.e., that the orbitals are biorthogonal. (A corresponding claim where the 
roles of and are reversed is equivalent.) 

Biorthogonality is equivalent to 

Wpv pnI^h,— ><in) = ^pi,<ii ' ' ' ^pn,in 
for the Slater determinant bases of V and V, given by 

\4>pr--p N ) — c pi c p2 ' ' ' c pn\~ ) anc ^ (^gi—gjvl = ( — \c-q N c<iN-i ' ' *5gu 

respectively. (We assume p± < p-i < ■ • ■ and q\ < q<i < ■ ■ ■ .) 

To prove the claim, suppose that (0 p >\<f p ') = 0. Then clearly, if = \4> P ', P2 --- ,p N ) , no 
£ V is non-orthogonal to Conversely, suppose that no such p' exists. Then, given 
an arbitrary 7^ |W) £ V, (<p pu ... , PN \^) must be nonzero for some Pi, - •■ ,Pn- This proves 
the claim. 

The creation operators are defined using field creation and annihilation operators as 

4-/*K. W «)td. (8a) 

and 

c p = (p p (x)i/)(x)dx. (8b) 
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The biorthogonality condition implies an anticommutator relation 

{c p , c\) = c p c\ + cjcp = {(p p \if q ) = S pq , (9) 

proven by inserting the definitions f!5a|) and f )8b|) . Thus, Wick's theorem {9} [ll| holds in its 
usual form, simply replacing c p with the operator c p . 

In standard CC theory, and virtually every other manybody method, f p {x)* = <f p {x), so 
that V = V^. We stress that the relaxation of this requirement allows for a complex analytic 
functional, which is essential for the bivariational principle. 



B. Excitation operators 

The orbitals are divided into occupied (the N first) and virtual orbitals (the L — N last). 
By common convention, indices i, j, k etc. denote occupied orbitals, while indices a, b, c, 
etc. denote virtual orbitals. The terminology comes from the fact that CC can be considered 
a perturbational scheme, where the exact wavefunction is written as 

l^) = e T |0), (10) 

where \4>) is a reference zeroth order approximation Slater determinant 

|0) = c\c\---c ] N \-), 

and T is an operator on the form 

t = r ^ + a* E T H 5 4tj + ■■■■ 

ai abij 

The operator c\ci destroys a particle in an occupied orbital (creates a hole) and creates a 
"virtual" particle above the Fermi sea defined by \<f>). The first sum on the right-hand-side 
is a singles excitation operator, while the second sum is a doubles excitation operator, and 
so on. A general n-fold excitation operator is on the form 

T = — r «i--:«n r t p. ... r t 

It is easy to see that without loss of generality the amplitudes can be taken to be anti- 
symmetric, which is the reason for the combinatorial prefactor. Importantly, since the occu- 
pied and virtual orbitals are disjoint sets, all excitation operators commute. It is convenient 
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to introduce a generic index /x for the excitations, i.e., 

where X M is a shorthand for X l a = c^Ci, X l ^ b = X l a X J b , and so on. Note that in the latter 
expansion, only linearly independent excitations are included, eliminating the combinatorial 
factors. 

It is worthwhile to note, that as an operator, T depends on both the amplitudes r = (r M ) 
and the orbitals through the X^, that is, T = T(r, It standard CC theory, one usually 

thinks of T as the primary unknown, since the orbitals are fixed, and since dependence on 
t is linear and one-to-one. In the OATDCC theory we must be careful, for example when 
computing dT/dt. Usually, however, there should be no danger of confusion when, for 
brevity, we suppress the parameter dependence of excitation operators. 

We also note, that even though X^ depends explicitly on the dual orbitals through the 
appearance of q, the function = X^tfi) does not: c\ is only responsible for removing (fi 
(not <pi\) from a determinant. The Slater determinants are easily seen to form a basis 
for V together with 

We define de-excitation operators as operators on the form 

S = J2 ^ = Yl + 212 E a ab4cac]c b + ■ ■ ■ . 

fi ai abij 

The term "de-excitation operator" has an obvious interpretation, but note however that 
these operators excite bra determinants. In particular, 

(0"|<M = $\X»X v \<f>) = 5*. 

Correspondingly, excitation operators de-excite bra states. 

We now comment on the form of the Hamiltonian in second quantization. For simplicity, 
we assume that the Hamiltonian contains at most two-body forces, i.e., 

N 1 

in first quantization, where h(i) is an operator acting only on the degrees of freedom for 
particle i, and u(i,j) acts only on the degrees of freedom of the pair For molecular 

electronic systems in the Born-Oppenheimer approximation, h(i) is the sum of kinetic energy 

12 



and the nuclear attraction potential, while u(i,j) is the Coulomb repulsion between electrons 
i and j. Suppose now II is the projection operator 

n = |0)<0| + ^|<?v)<n 

which acts as identity on V[<§] and V[$]: For any G V, II|^) = |$), and for any (^'| G V, 
(\l/'|n = However, II 7^ 11^ so it is not an orthogonal projector. We now have 

($'\H\y) = (tf'|mra|tf), 

where 

UHU = J2(<p p \h\<Pq)clc g + - ^2(<p p <p r \u\<p q <p s ) AS clclc s C q 
pq prqs 

pq prqs 

The two-body integrals are anti-symmetrized according to the standard in CC theory and 
are given by 

(ip p(p r \u\(f q (p s ) as = {(p p (p r \u\ip q ip s ) - ((p p (p r \u\(p s ip q ), (12) 
(0 p if r \u\if q if s ) = J <p p (x)<p r (y)u(x,y)<p q (x)<p s (y) dxdy (13) 

It is important to note, that unless the orbitals are complete, ILHTI 7^ H. 

Similar considerations as the above also hold for arbitrary one- and two-body operators. 

C. From variational to bivariational CC 

Having discussed orbitals and operator expressions using second quantization, we now 
turn to the CC ansatz. It is a fundamental fact of CC theory that any wavefunction |$) 6 V 
normalized according to (<ft\^) = 1 can be written on the form fflUj) . i.e., the exponential 
ansatz is covers the whole of the discrete Hilbert space V. To see this, we simply observe 
that exp(T) = 1 + A, where A is a new excitation operator. By writing T = T\ + T2 + • • • 
and A = Ai + A2 + ■ ■ ■ , A and T can be compared term-by term, giving explicit formulae for 
Tfc in terms of Ag, i < k, and vice versa. Since any \^) G V with (<f)\^) = 1 can be written 
I*) = \4>) + A\$), the result follows. 
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Similarly, we have (^'\ = {4>\e T ' for any (^'| normalized according to (&'\(j>) = 1, where 

A" 

is a de-excitation operator. Note that the exponential parametrization results hold for any 
choice of biorthogonal orbitals. 

We note that a truncation in A instead of T at, say, A = A\ + A2 gives a linear 
parametrization which defines the CI singles and doubles ansatz, CISD. [Compare also with 
Eqn. ©.] 

The bivariational expectation Sjj now reads 

£hKt ^^M, (14) 

where we note that the functional dependence on the orbitals is implicit in the reference 
determinants and T and T' . 

We now make the observation, that (oj\4>) = 1, where 

U = (aie^e^Y 1 (d>\e T 'e T , 



implying that there exists a de-excitation operator S = Ylu a ^^ sucn that (^1 = (^l^ '• 
From this we obtain 

(*'| = (0|e T 'e T |0)(w|e- T = (0|e T, e T |0)(0|eV T = (0|e T 'e T |0)(§|. 

Inserting this into ( !T4|) . we get rid of the denominator, viz, 

S H [(T, r, $, $] = (V\H\V) = (4>\e s e- T He T \<f>). (15) 

We stress that there is no loss of generalization in these manipulations. 

We perform a further change of variables. Writing e s = I + S + ■ ■ ■ = / + A with 
A = partially transforming from exponential to linear parametrization of 

£ h [\,t,$,$] = (0|(/ + A)e- T #e r |0>. (16) 

Disregarding the dependence on $ and $, this functional is the well-known CC expectation 
functional jlo|, I2I 13 j - The usual interpretation for A M is as Lagrange multipliers for a 



constrained minimization of the energy which is equivalent to the standard CC equations 
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15] . In our case, however, it is interpreted as part of the parametrization of the approximate 
dual wavefunction which enters the expectation functional. should be treated on equal 
footing with r M : they are equally important. 

All three functionals (j!4p . (TT5"j) . and (FTO]) are equivalent to the multi-configurational 
Hartree-Fock (MCHF) functional. The fundamental approximation idea in CC theory is 
now to truncate the expansion for T and A (or S if the functional ( Tl5l) is used) at a finite 
excitation level. For example, for the coupled-cluster singles and doubles (CCSD) approxi- 
mation, one takes 

T^T 1 + T 2 and A » Ai + A 2 , 

neglecting amplitudes with n > 2. In the following discussion, the truncation level n should 
be considered a parameter of the ansatz. Under such a truncation, the three CC functionals 
Eh[t' , r, $, $], £h[o~, t, $, $] and £f/[A, r, $, $] are no longer equivalent. 

Suppose for the moment that the orbitals are held fixed and orthonormal. What methods 
do the three functionals ( IT^jl - fTloT) define? The functional ( TH1) corresponds to a variational 
coupled-cluster theory (VCC) within the chosen basis and is abandoned for reasons that will 
soon be apparent. The functional (Tl5|) defines the so-called extended CC (ECC) functional 



10|, and defines an alternative approach to the standard CC defined by the CC Lagrangian 
( fT6|) . In fact, it may be viewed as a more "natural" version of CC since is parametrized 
exponentially, i.e., size-consistently. 

In the remainder of the paper, we focus on the standard CC Lagrangian ( fl6l) . but with the 
orbitals are included as variational parameters. We correspondingly refer to the functional 
as the "OATDCC expectation functional" in the rest of the paper. 

For the reader not acquainted with CC theory, it may be hard to see that we have 
actually simplified matters with the CC Lagrangian or ECC functional compared to simply 
considering the VCC expectation value functional, with T' — T\ which would be the default 
approach of optimizing the energy using any ansatz. The problem is that while the ansatz 
scales polynomially, the evaluation of (0|e Tt ife T |0)/ (0|e Tt e T |0) does not. It also seems hard 



to find size-extensive approximations of finite order in T [32]. 

On the other hand, one of the basic observations of CC is that Baker-Campbell-Hausdorff 
(BCH) expansion of the similarity transform exp(— T)H exp(T) truncates identically after 
a finite number of terms, regardless of the number of particles present. In fact, for a Hamil- 

15 



tonian with at most two-body potentials only terms up to fourth order are nonzero, 

4 i 



n=l 



where [H, T] n = [[H, T] n _i, T] is the n-fold nested commutator. This no less than remarkable 
fact can be seen from [[c|, X^], X v ] = [[c p , X^], X„] = (verified by direct computation) and 
by the fact that H is a fourth order polynomial in the creation and annihilation operators. 

It follows that £#[A, T, $, $] is a fourth order polynomial in r = (r M ) and linear in 
A = (A^). This polynomial can be evaluated using Wick's theorem. As this involves an 
agonizing amount of algebra, an alternative approach resides in the use of Feynman graphs 



to simplify Wick's theorem 0, [l5|, or in the use of computer algebra software 
latter approach has become more common in the last years and is utilized here, see |A] 



33|. The 



D. The OATDCC action functional 

Having established the form of the OATDCC expectation functional, we now turn to the 
evaluation of the corresponding action-like functional S defining the Schrodinger dynamics, 

S[X, r, $, $] = £(4>\(1 + A)e" T (i^ - # )e T |0> dt. 

To evaluate the explicit functional dependence on the time derivatives of r and $, we must 
compute J^|\l/) = J^e T |0). To this end, we use the expansion 

|tf> = |0) + J>1^>. A" = A"(r) = (r\e T \<f>), 

where the summation is not truncated at any level. Since Wick's theorem only uses the 
anti-commutator fl9]), the coefficients = A^ij) do not depend explicitly on the orbitals, 
only on the (possibly truncated) amplitudes r. Moreover, we compute the derivative of a 
Slater determinant via 

d 

— c' c> • • • c' I— ) = • ■ • c' I— ) 4- c' • • - c* I— ) + 

= {^-< ^<f^j C Pl C P2 " ' C PN I ~~ ) • 



1(3 



This implies 



|i*> = E (s"" (T) ) w + (e^J w + E^m (E4 5 «J 

e^+(e* 



(17) 



In the last step we used 
For the time-derivative part of the functional integrand, we now get 



iM0|(l + A)e- T |e r | 



ih{4>\{l + K^)e~ T \Y,t v X v + Dj e 1 

ih^x^ + ih(4>\(i + A)e^ T n J on e T |0). 



The projected operator ILDII is easily computed using 

\(p p ) = {P + Q)\^ p ) = ^ \<Pq){<Pq\<fp) + Q\<Pp) 

q 

where P = is the (oblique) projector onto single-particle space, and Q = 1 — P. It 
follows that 

ILDII = ^{<p p \<p q )clc q = D . 



pq 



Finally, we obtain 



S[A,r,$,$] = f ifiy]A/-5 ff _ ifiDo [A 1 T I $,$]df 
Jo 



T 



where 



ih\,f» + pl{h\ - ihrf q ) + dt 



pl = pl(X,T)^m + A)e- T clc q e T \ ( f>), 
Pit = PU\ t) = (01(1 + k)e- T c\c\~c s c q e T \ 
K = h pq {®,$) = (0 p \h\<p q ), 



(18a) 
(18b) 
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and 



<s = = (<fp<fr\u\<fi q <fi s ) 



AS- 



In Eqn. (1 1 8b j) we introduced the Einstein summation convention over repeated indices of 
opposite vertical placement. This greatly simplifies the algebraic manipulation of CC ex- 
pressions, see \M (Strictly speaking, we then should introduce a similar index placement on 
the orbitals, e.g., <p p instead of (p p , and cP instead of c p . However, we find this a too great 
departure from the standard, so we, keep a somewhat inconsistent notation for simplicity.) 

The quantities p q p and p qs r (whose index placement should be noted) are the CC reduced 
one- and two-particle density matrices, respectively, and are not explicitly dependent on 
the orbitals, since they are evaluated using Wick's theorem depending only the fundamental 
anti-commutator. They therefore only depend on the amplitudes. Similarly, the one-particle 
integrals M, r\ v and the two-particle integrals u p Z. only depend on the orbitals. These facts 
dramatically simplify the application of the evaluation of 5S. 



E. Standard CC and linear response 

For clarity, we briefly discuss the standard CC method for the ground state problem 
and the time-dependent generalization, i.e., the equations of motion used in linear response 
theory jl^ . For standard CC, a fixed set of orthonormal orbitals are chosen, i.e., 

(0p\ = \Vp)\ {c P ,cl}=5 pq . 

The orbitals are usually but not necessarily chosen to be the Hartree-Fock orbitals for the 
given (bound) system. 

Since the orbitals are fixed, the only parameters in the expectation value functional are 
the amplitudes r and A, 

S[X,t] = (<P\(l + A)e- T He T \(j)) 

= E cc [t] + \Me- T He T \(P), E cc [r] = (0|#e T |0>, (19) 

where we have used (</>\T = 0. Computing the derivatives with respect to A M and equating 
to zero, we get the stationary conditions 

(<P^\e- T He T \<P) = 0, V^GZ, 

18 



where the set X contains all the amplitudes in the desired approximation, say CCSD. Note 
that we could arrive at this equation by simply similarity-transforming the Schrodinger 
equation, 

e- T He T \4>) = E\0), 

and projecting against 

By projection onto (0| we get 

E = (<l>\He T \<i>) = Ecc[t], 

which is of course exact within V = V if the amplitudes are not truncated. Truncations give 
an approximate similarity transformed Schrodinger equation, and the interpretation of A M 
as Lagrange multipliers for the constrained minimization of Ecc is then apparent. 

In the CC literature, the coefficients A M were originally introduced jn order to compute 



expectation values consistent with the Hellmann-Feynman theorem 
one seeks an expectation value functional defined by the criterion 



11 
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i.e., that 



(A) = §-E CC (e] 



e=0 



where i?cc( e ) is the energy eigenvalue approximation found with the CC equations for the 
perturbed Hamiltonian H + eA. (Recall that (A) is the first order perturbation in Rayleigh- 
Schrodinger perturbation theory.) It turns out, that (A) = £a[A,t], where (A, r) is the 
critical point of the functional (TT9]) . i.e., the unperturbed CC solution. 

For the time-dependent case, one simply similarity transforms the time- dependent 
Schrodinger equation, obtaining 

iht\<p) = e- T He T \<P), 

and via projection, 

iftH* = (<f )ll \e~ T He T \ ( f ) }. 

This is exact (within the space V) for the untruncated ansatz, but motivates the use of this 
equation in the truncated case as well. In that case, it is easily obtained from the stationary 
conditions of the CC Lagrangian 

S[\,t}= [ (0|(l + A)e- T (ift|--#)e T |0) dt = [ ih\^ -S[\,r] dt. 
Jo Jo 
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These are 







S[X,t] = (^\e- T He T \<P) 



(20a) 







£[X,r] = (4>\(l + A)e~ T [H,X,}e T \ ( p) 



(20b) 



as is easily verified; seeO Expectation values are computed using £U[A, t] as in the stationary- 
case. 

IV. THE OATDCC EQUATIONS OF MOTION 
A. Parametric redundancy 

In order to derive the equations of motion for r = (r M ), A = (A M ), $ = (tp p ) and 
$ = (if p ) given by the stationary condition 5S = 0, we must, in the functional S, insert 
all independent variations of the parameters. However, for a given CC wavefunction pair 
u = ((^l, I 1 !/)) G A4, there are many choices of the amplitudes and orbitals that give 
the same pair. Correspondingly, not all variations give independent equations, and this 
ambiguity must be eliminated. 

For example, it is well-known in CC theory that occupied and virtual orbitals may be 
rotated among themselves with a corresponding inverse rotation of the amplitudes, such 
that the wavefunction is actually invariant to within a constant factor. 

To formalize this somewhat, we write the collected parameters as a point z in a manifold 
Af, i.e., z = (A, r, $, <3>) G Af. The CC ansatz is then given by a many-to-one mapping f(z), 



where we for the action functional explicitly have written out the dependence on the history 
for clarity. Suppose we identify a Lie group G such that 



M, 



inducing functionals 



E[z]=£ H [f(z)}, and S[z(-)} = S[f(z(-))}, 



E[G oz}= E[z] for all G G G, z G M 
S[G(-) o z(-)} = S[z(-)} + C for all G(t) G G, z(t) G Af 



(21b) 




20 



N = \z= U,T,<t>,4>) 



orbit of G 



z + 6z' 




f(z) 




FIG. 1. Illustration of the parameter manifold M and the CC wavefunction manifold Ai. An 
element u £ At is given as the image u = f(z), albeit in a non-unique way. The pre-image f (x) 
is the orbit of a Lie group G acting on some z, illustrated as a closed curves on passing through z. 
Thus, there is a manifold of displacements bz giving rise to the same displacement 5u. Specifying 
a gauge selects a unique 5z (a unique G) for each Su. 

with C being a constant (only dependent on G). Equations fl2TT) states that the CC function- 
al are invariant under the action of the Lie group. Let z(t) be given. We now observe, that 
there is a continuum of histories G(t) °z(t), one for every choice of G(-), that corresponds to 
the same integrand. So any change 5z(t) = G(t)z(t) with G(t) infinitesimal will not change 
the value of the action functional. Such variations of z(t) will lead to redundant equations. 
We comment that infinitesimal G(t) is on the form 1 + g(t), where g(t) is in the Lie algebra 



From the invariance of the functionals, the physics predicted by G o z is the same as that 
predicted by z. The physical solution is unique, but the parameters are not. This situation 
is similar to the one in gauge field theories, so the elimination of these extra degrees of 
freedom in the parameters is called a gauge choice. In Figure [1] the relationship between z 
and u = f(z) is illustrated, along with the concept of the invariance under the action of G. 

The group G describing mixing of occupied and virtual orbitals separately consists of all 
block diagonal and invertible matrices on the form 



i.e, the matrix elements G ia = G ai = 0. The action on z is denned as follows: The orbitals 



of G. 
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are transformed according to 



$ — > $G, $ — ► G~ 1( l>, (22) 
which preserves biorthogonality. Equation 022 p is equivalent to the transformation 

C p ^ C lGqpi Cp ^ G pq Cg. 

1 Q 

This gives the following transformation on X % a : 

We define the transformation of the amplitudes by 

— ►£ G W G a'i (23a) 
i'o' 

T ii ^ E G H lG n lT i'j' Ga^cPb'b ( 23b ) 

i'j'a'b' 

K^Y, G ^' G ^ ( 23c ) 

i'a' 

*L — ► E GoafG w r%G^Gji (23d) 

with corresponding expressions for higher order excitations. This completes the definition 
of the action G o z. 

As a consequence of the transformation, T — >■ T, A — y A, while |</>) — y g\<p) and 
(0| — y (</>|g _1 , where g = det(Cry). Clearly, E[z] = E[G o z]. As for the time-dependent 
functional S'[z(-)] it remains to check the part containing the time-derivative, i.e., 

f T ~ r) f T ~ 8 f T ~ r) <9 In n 

I <*!«!*> At ^L ^"V 1 *^/ <*^*> + ^r di - 

so the integrand gains only a total time derivative, 

S[z{-)\ — ► S[G{-) o *(•)] = S[z(-)\ + ihHg(T)/g(0)]. 

The last term is a constant with respect to variations vanishing at the end points of < t < 
T. 

We now proceed to choose a gauge for the orbital rotation group. Suppose we have 
a solution candidate Zo(t), i.e., 5S[zo(-)] = f° r an y variation of z (t). For any choice 
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G(t) G G, z[t) = G(t) o zo(t) is also a solution. To fix a unique solution z(t), we need to 
find a differential equation for G(t) in terms of zo(t) with a unique solution. To this end, 
consider the action of G(t) on the orbitals, i.e., = %(t)G(t) and = G(t)~ 1( lo(t). 
Writing rj (t) = & (t)& (t) an d = w e have 

6?(^(<)=%(<)G(0 + G(t). 

We are free to choose the matrix elements G(t)ij and G(f) a b of the nonzero blocks G(t) occ 
and G(t) v i T , respectively, of G(t) at will. Let g(t)ij and g(t) a b be the matrix elements of some 
arbitrary matrices g(t) occ and g(t) v i r , respectively, and require 

G(t) occ = G(t) occ g 

(t) occ occ 

G(t) vir = G(t) vir g(t)vi T - r]o(t) vh -G(t) viT 

which defines G(t) uniquely. Then, 

v(t)occ = g(t)occ, and /?(t) vir = g(t) vi r 

so that any choice of g(t) occ and g(t) v i r , i.e., choice of gauge, implies a specific choice of 
transformation G(t), and hence a solution representant z(t), all being equivalent 35]. 
The simplest choice is probably g(t) occ = and g(t) vir = 0, such that 





17(0 = $(t)$(t) 



iia\ 





To derive the differential equations corresponding to this gauge, one needs to perform all pos- 
sible variations Sz(t) adhering to the constraint (<M$) occ / vir = 0, i.e., (<fii\Sipj) = (<p a \5ipb) = 
0. 

Note that G(t) is simply a theoretical device that describes the continuum of solutions 
z(t). Using this device we derived conditions on z(t) and the variations 5z(t) that picks 
exactly one of these solutions. G(t) will not actually be solved for, since it has no physical 
value. 

However, G as defined above is not the largest group leaving the functionals invariant. It 
turns out that the singles amplitudes rf may be completely eliminated as well, corresponding 
to orbital transformations on the form 



G 
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To see this, note that 
l*> = e T k 



e T - Tl e Tl l 



e Tl c j N e~ Tl 



)!->. 



and that the transformation cl — > exp(Ti)d,exp(— Ti) is equivalent to (see IC]) 







(o 0^ 






r 




$ - 








->■ e 







The dual state (W| is invariant: 

= (0|(1 + A)e" T — »• (0|(l + e Tl Ae^ Tl )e- T+Tl 



<*l, 



where we have used that (</>|e Tl = The bivariational functionals are then invariant under 
the described transformation, and we may set rf = 0. 

Observe, that we did not transform the amplitudes X l a . It is tempting to assume that a 
transformation on the form 

^0 at} 



G = e 



may achieve an elimination of A^. This is not the case. There is no transformation 
with the same truncation level that compensates for the transformation 

m = e T \<p) 



e Sl e T 'e~ Sl e Sl \ 



e Sl e T '| 



i.e., G induces higher-order excitations in T. The bivariational functionals therefore cannot 
be made invariant under the transformation, so A^ cannot be eliminated in this way. 

However, it is easy to see, that if we do include A^ as parameters the equations of motion 
will be overdetermined. The presence of T\ is compensated by the freely varying orbitals, 
but the same is not true for A x . We conclude that in the orbital- adaptive CC, the dual 
state (\&| would have more parameters than |\&). Correspondingly, (\I/| and |W) would not 
be one-to-one. 

From now on, we therefore set both rf and A^ identically equal to zero, in effect starting 
with a coupled-cluster doubles, triples, etc, model with adaptive orbitals. 



24 



B. Derivation of equations of motion 



Performing the variations with respect to the amplitudes is straightforward (see IC]) . and 
leads to the equations 

^ = ^H-ihD [X: r, $, $] = (^\e- T (H - ihD )e T \ ( f ) } (24a) 

-ihX^ = ■£ jI £ H - iRDo [\,T,$,$] = (0|(l + A)e- T [^-i^ o ,^]e T |0) (24b) 

which must hold for all fi G X included in the approximation. Equations ( l24"j) are identical 
to Equations f!2Up except for the presence of D Q due to the changing orbitals. 

Performing the variations with respect to the orbitals, we observe that an arbitrary 
variation of the one-body function \(p p ) = \(p p (t)) (where we use ket notation for clarity) can 
be written 

S\<p p ) = PS\ip p ) + Q5\(p p ), 

with P = = J2 q \ip q ){<Pq\ and Q — 1 — P. Due to the gauge choice rf- = rj^ = 0, we have 
for the occupied and virtual orbitals 

8\<Pi) = y £ i 4\<Pb) + Q8\<Pi) 

b 

and 

j 

respectively, where ef are arbitrary independent functions of time, and Q5\<pt) is completely 
arbitrary and independent from the e". We may therefore perform the variations in two 
stages: First, we set 6\<fii) = e(t)\tp a ) for each a,i separately (then exchange % and a), and 
then finally we set 5\(p p ) = \6) = Q\6). 

Beginning with the P-part of the variations, we note that the variations in (<p p \ are linked 
to those of \(f p ) due to the biorthogonality constraint. We get 

= 6(0 a \<fi) = {6<p a \<Pi) + e(t) <^> 6(<fi a \ = -e(t){(fii\. 

It is useful to consider an arbitrary \u) and (v\, whose variations are 

5\u) = ec^ a Ci\u) and 6{v\ = — e(v\c\ci, 
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so that we have the general equation 

5(u\A\v) = e(u\[A,clci\\v) + (u\(5A)\v), 

where A is any operator that may depend explicitly on the orbitals. (If A does not depend 
on the orbitals, then 5 A = 0.) For example, the Hamiltonian does not depend explicitly on 
the orbitals, while 

5 D = 5 c% = i ic l + e 4) d * ~ ed l d i = kc %- 
<? 

We also note that 



<^f^|X M |*) = 0, 



since the expectation value of X M only depends on the amplitudes and not the orbitals. We 
compute the variation in S: 



5S = df (ylihDo + ihJ^^X^- H\V) dt 

= / e(t)^\[ihD -H,clcm + ^Ht)^\cid^)dt (25) 
Jo 

= J e(f) ((¥|[iftD - H,clcm - ihpl) dt. 

Requiring 5S = for all e(t) implies that the integrand must vanish. Using the gauge 
condition, D becomes 

jb jb 

and only the first sum survives in the commutator in Eqn. (1251) . We obtain the equation 



bj 

which is a linear equation for rf- = (<p b \<pj). It is readily verified that 

Pa = (*\<im = (01(1 + A)e- T Xfe T |0) = A* = 0, 



and that 



[ c p^<?> 4 C i] ^a C p C i 3pC a Cq (27a) 
CsCqj C^ a Ci\ S^cLc^CsCi 8 a C!pC! r CqCi -\- 5 r C^ a CpC s Cq 5pC^ a C^,C s Cq. (27b) 
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The coefficient matrix of Eqn. ( |26|) becomes 

In total, we get a linear equation for rf- that reads 

1 r £4 



(28) 



fej p (jr L P rs 

where we used Eqns. (l2~?]) and the anti-symmetry 



s v pr 
pr "'as 



Far u qs 



rqs 



-U 



pr 
sq ■ 



Y,P a p r<-J2p 



ir qs 



prs 



rqs 



We now turn to the variation 5\<p a ) = e(t)\ipi), which implies 5(<pi\ = —e(t)((p a \. The 
calculation is completely analogous to the previous case, so we simply state the result, 

bj p q 

Unlike p l a , the coefficients pi do not vanish identically, except for in the doubles only ap- 
proximation to be considered in Section IIV CI 

Having derived the equations of motion for the P-part of ip p (and therefore also (p p ), 
we now turn to the Q-part. We perform an arbitrary variation 5({p p >\ = (6\ = (9\Q, i.e., 
(0\(p q ) = for all q. These variations are therefore independent from the corresponding 
variations in \<p q ). 

It is convenient to use S on the form (|18bl) to derive the variational equations. Recall 
that the reduced density matrix elements p q v and p qs r are independent of the orbitals, since 
they are computed only using the anti-commutator (Q. Thus, the only quantities that vary 
are the one- and two-body matrix elements h p , r] P and u p q r s : 



5S = 6 



E^K-^)-iE« dt 

pq prqs 

[ J^Pl'i 9 ^^ - h )M - \^P¥r( G( Pr\u\ip q ips) dt 

qrs 

pl,(ih\cpq) - h\<fg)) - E^p'r( • <fr\u\<fq<fs) 



(29) 



dt 



q qrs 

From line 1 to line 2 we used the symmetries u pr s = u rp and p qs r = p^, and the two-electron 
integrals in the last two lines are not anti-symmetrized. We define mean-field potentials WJ 
by 

( • <p r \u\<p q tp s ) = / <p r (x')u(x,x')<p q (x)<p s (x') dx = WZ\(p q ). 
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Since Eqn. ( 129|) must hold for all (9\ = (9\Q, we get the equation 

d 



Performing the variation 5|<^ 9 ) = \9) = Q\9) is completely analogous, and gives 

-' ih Y,pi (J^i) Q = E^>^ + Y,p q ^p\ w :Q for a11 9> 

p ^ ' p prs 

where the minus sign comes from integration by parts. 



C. The doubles approximation: OATDCCD 

The simplest non-trivial OATDCC case is the doubles approximation (OATDCCD). The 
wavefunction parameters are, in addition to the orbitals $ and $, the amplitudes r = (r? fe ) 
and A = (A^)- In [A] a complete listing of the algebraic expressions needed to evaluate 
the equations of motion is given. In particular, the only nonzero elements of the one- 
body reduced density matrix are p{ and p b a , which simplifies the equations of motion. The 
amplitude equations read 



ihf r 







ab 

~ [hX ab = o-^fe 



^[A,T,<M] 



l%\e- T He T > 



S H [X,r,<!>^} = mi + A)e- T [H,XZ]e^ 



The P-space orbital equations read 

bj j b 

bj b j 

Finally, the Q _s P ace orbital equations are 





prs 

E 

prs 



rpr as 






rqs 


o as u pr 

rpr "'ts 






rqs 



ar qs 



qs ar 
ir qs 



q q qrs 

ih Y.pi (J^i) Q = E^w +T,p q pr&p\wiQ. 



(30a) 
(30b) 



(30c) 
(30d) 



(30e) 
(30f) 
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The coefficients A 1 ^ are defined in Eqn. (128]) . The right-hand sides of Equations fl30al) and 
(130bl) . which are polynomials in terms of rff, A^ 6 , h p q and u p q r s , can be found in the Appendix. 
Note that the right-hand sides of the equations are identical to the one used in standard 
CCD calculations for the ground state energy, since the operator D is eliminated due to 
Pi=Pa = 0. Existing computer codes may be helpful for implementations. 

Since D drops from (I30aj) and fl30bl) . the right hand sides can be evaluated independently 
of Eqns. (I30c[) to (130f|) . Note that to evaluate $ and r] must be solved for in addition to 
Q\<p q ) and (ip p \Q. The results are assembled according to 

\<Pq) = {P + Q)\Vq) = M^P^l) +Q\Vq) = ^ V^lfp) + Q\<Pq) 

P P 

(<Pp\ = (<P P \(P + Q) = ^2{<Pp\<Pq){<Pq\ + (<Pq\Q = ~ ^ < (<A? I + (<Pp\Q 

q q 

We note that the Q-space equation for \<p p ) is formally identical to the orbital equation of 
MCTDHF (see[B]), and the equation for (<p p \ is formally identical to the complex conjugate. 
However, in the CC case the matrices u^., h v q and pp" are not exactly Hermitian. Therefore, 
the two equations are only complex conjugates of each other to within an approximation, 
and both must be propagated. 

We will now consider the computational cost of evaluating the time derivatives in a 
computer implementation. We will in the following assume a grid-based discretization of 
single-particle space using in total iV& points. In particular, integrals are evaluated as sums 
with Nb elements. 

We consider first the computation of the amplitude equations f|30ap and f!30bj) . assuming 



that h v etc are available in computer memory. It is easy to see, from Equations (1A70 and 
( 1A12|) . that by brute-force summation the worst-scaling terms require 0(N (L — N) 4 ) = 
0(L 8 ) operations for computing the totality of derivatives, which is a conservative estimate 
since N < L. Existing CC codes typically reduce this to 0(L 6 ) by clever use of intermediate 



variables 15]. 



The P-space orbital equations fl30cl) and f!30dj) are linear equations where a vector of 
dimension 0(L 2 ) is to be solved for. This requires at most 0(L 6 ) operations. The right- 
hand side is dominated by the two-body terms, which cost 0(L 5 ) in total to compute. 

We next turn to the Q-space orbital equations (I30el) and (130fl) . which can be viewed as 
differential equations for matrices of dimension JV& x L and L x Nt, respectively. The cost 
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analysis is identical for the two. The matrix p q needs to be inverted, a step of at most 0(L 3 ) 
cost. Multiplying Eqn. ( 13 Pep by the inverse matrix elements of p q v shows that the Q-part 
of ih\(p p ) can be computed as the sum of Qh\tp p ) and a two-body mean- field term which 
clearly dominates the computation. The cost of this term is 0(L 3 Nb) plus 0(L 2 Nb) for the 
multiplication of the result with p _1 . 

Unlike standard CC calculations, the one- and two-electron integrals hP and u pr s must be 
updated at each time t. This is similar to the situation in MCTDHF theory. Moreover, 
equations (I30cl) - fl30fl) are formulated in terms of the reduced one- and two-electron matrices 
which need to be computed. 

The matrices p q and p qr s cost less than the evaluation of the orbital equations right-hand 
sides in total, and h p is relatively cheap to compute. However, the two-electron integrals and 
the mean-field functions W£ are costly. The computation of all the mean-fields, which are 
local functions, costs 0(L 2 N 2 ). Since u pq s = (0 P \ W^\ip q ), the computation of the two-electron 
integrals costs an additional 0(L 4 Nb) operations. 

The computation of W r s is in fact very expensive, being similar in cost to computing 
two-particle integrals, i.e., six-dimensional integrals in realistic calculations. This problem 
is ubiquitous for all time-dependent mean-field calculations, and a common approach is to 
employ some low-rank expansion for the interaction potential u(x,x'), which needs to be 
sampled at the grid points k = 1, • • • , N b . The resulting matrix v(£,k,£,k') is symmetric, 
with eigenvalue decomposition 

N b 

V(£k,£k') = ^2 X mfm(Ck)fm(Ck'), (31) 
m=l 

where f m is the eigenvector belonging to A m , the latter arranged in decreasing order. The 
optimal (in the 2-norm) M-term approximation to is then obtained by truncating 

Eqn. (13TT) after M terms. Oftentimes, only a small number M Nf, terms are needed. 
Moreover, if is increased, M may typically be held fixed. The mean-fields then become 

M 

{£,k){<fir\v m \ip s ), 

m=l 

which reduces the cost of computing W£ to 0(NbL 2 ) for the inner products plus 0(L 2 NbM) 
for the summation, reducing the cost proportionally to the fraction M/Nf, of modes included 
in Eqn. fl3TJ). 
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To sum up, we see that the cost of evaluating the right-hand sides of the equations of 
motion is dominated by the computation of the WJ and u^ r s and of the evaluation of the 
amplitude equations, costing 0(N 4 (L — iV) 4 ) (if no optimization is done), and being the 
only terms that increase in complexity with the number N of particles. This should be 
contrasted to MCTDHF calculations, where the number of amplitudes grow exponentially 
with N. 

V. FURTHER PROPERTIES OF OATDCC 
A. Relations to other methods 

It is instructive to consider special cases of the OATDCC method and relate these to 
other, well-known wavefunction approximations. 

In the case where all excitation levels are included, we have seen that the CC ansatz 
becomes the FCI ansatz within the chosen basis. Therefore, the MCTDHF and OATDCC 
methods are equivalent in this limit. In [B] an explicit derivation of MCTDHF using the 
bivariational principle is carried out. Note, however, that the gauge conditions are different, 
i.e., the orbitals produced are not identical. This stems from the fact that the CC wave- 
function is normalized according to = 1 at all times. However, the orbitals in the two 
methods span the same space, and the wavefunctions are identical except for normalization. 
At each time t, the OATDCC and MCTDHF wavefunctions produce the same expectation 
value functional for any observable. 

At the other end of the hierarchy, we find the trivial case where there are no amplitudes 
at all, i.e., we take L = N and A = T = 0. In that case, the OATDCC energy expectation 
functional becomes 



£„[*,*] = WW) = M 



\ pq pqrs 



uP qs C p C r^Cq \<f>) 




V 



pr 



This is the Hartree-Fock energy functional (when $ = $ H ). Thus, the conditions for 58 h = 
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are the Hartree-Fock equation and its complex conjugate. The action functional becomes 

= / (4>\{ihD -H)\<f>)dt 
Jo 

^ p p pr 

This is the time- dependent Hartree-Fock (TDHF) functional. Computing the variation with 
respect to (p q we get the TDHF equations of motion, and the variation with respect to (p q 
gives the complex conjugate, showing that the T = A = case is indeed equivalent to 
TDHF. 

We also note that the OATDCCD approximation is equivalent to MCTDHF whenever 
N = 2. Moreover, some combinations of L and N also give equivalence, for example L = 
N + 2 since there are no triple excitations defined. 

Finally, in the absence of interactions, the Hamiltonian is a pure one-body Hamiltonian. 
One can easily show that the choice iM> = — iM> = and A M = f M = gives 

S[X, t, <3>, <3>] = and a stationary S. This is the exact solution to the dynamics for any 
initial condition. (The gauge conditions in this case are chosen differently from earlier: 

= (<Pi\<Pj) = (ft|# (1) bj)M and V a b = (<p a \HW\<p b )/iK.) 

In OATDCC the evolution of the orbitals are chosen to variationally optimize the action 
functional. As early as 1978, Hoodbhoy and Negele [19J discussed a time-dependent CC 
approach using an explicit dependence of time in the orthonormal single-particle functions. 
This would correspond to using the following functional to define the evolution: 

Sh-n[t, A] = / ih\^ - (0|(1 + A)e~ T (H - ihD (t))e q » dt. 
Jo 

The Dq operator simply is a correction in the standard CC Lagrangian due to a moving basis. 
Hoodbhoy and Negele suggested computing the time-dependence of $ using, say, TDHF, 
i.e., a TDHF calculation is first performed, and the output is fed into <Sh-n- However, 
this approach would have inferior approximation properties compared to OATDCC while 
at the same time being only marginally easier to evolve in time. To see this, consider the 
fact that the TDHF solution would depend only on the initial choice of the orbitals, and 
not on the full state at time t as in the OATDCC approach. The OATDCC orbitals will 
generally differ substantially from the TDHF solution, since their motion is computed from 
the wavefunctions at time t. The Hoodbhoy-Negele TDHF approach neglects the correlation 

32 



effects built into the wavefunction during the evolution. As for the computational cost, note 
that TDHF needs the computation of the two-particle integrals, so one gains very little, if 
nothing, at simplifying to TDHF for the orbitals. 



B. Approximation properties 

We now ask: what kinds of systems can we expect to be able to treat with OATDCC, 
and what systems cannot be expected to give good results? 

The usual CC ansatz (with fixed orbitals) is based on a single reference determinant, 
incorporating correlations through the cluster operator. For good results, the reference 
determinant should be a "large" part of the wavefunction. In the language of computational 
chemistry, dynamic correlation (which by definition is due to the interparticle interactions) 
must be dominating, while static correlation (arising from degeneracies in the spectrum) 
should be small. 

For a dynamical calculation we must correspondingly require that, for all £ > 0, the wave- 
function is a single-reference type state. This is reasonable whenever the initial condition is 
of such type: intuitively, the Hamiltonian cannot generate static correlation since the only 
source of correlation from dynamics is the interparticle interaction. This is actually observed 
in the numerical experiment in Section I VI I 

However, standard CC is known to perform adequately even in the presence of static 
correlation, which gives reason to believe that the same holds true for OATDCC calculations. 
In any case, OATDCCD should be much better than a singles-and doubles truncation of 
MCTDHF due to size-consistency, even though the latter is a "true" multiconfigurational 
method, where all the basis determinants are independent. 

It is worthwhile to note, that for systems with spin, the OATDCC ansatz is not an 
eigenfunction of the total spin; only of the total spin projection along some preselected 
direction in space z. This is, however, not a problem in general if H commutes with total 
spin: all expectation values for spin-independent observables are the same as for the properly 
"spin-symmetrized" wavefunction. 
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C. Non-feasibility of imaginary time relaxation 

We have not yet discussed choices of initial conditions for OATDCC calculations. Typi- 
cally, one would like to start in the ground state of the system under consideration, or a state 
closely related to this. Indeed, the orbital-adaptive CC ansatz could in principle be used for 
ground-state calculations in the first place, just like MCTDHF actually is a time-dependent 
version of the multi-configuration Hartree-Fock (MCHF) for computing eigenvalues of H. 

The ground state is a critical point for Eh, and in MCTDHF theory this can be computed 
using imaginary time propagation, that is to say, formally replacing the time t with —is; 
so-called Wick rotation. Asymptotically, as s — > oo, a critical point of the variational 
energy is obtained. This is a quite robust procedure for variational approximations. For 
non-variational methods like coupled-cluster, the situation is, literally, more complex. 

To see this, consider again local coordinates z(t) G C n . The equations of motion for z(t) 
are analytic in both z and t. Thus, the energy is an analytic function of t, which is conserved, 
d£[z(t)]/dt = 0. Thus, the energy is a constant analytic function, even for complex t. Thus, 
the energy will not decay exponentially when the system is propagated in imaginary time. 

We conclude that computing the ground state using imaginary time propagation is not 
straightforward and requires separate study. Instead, quasi-Newton schemes like those al- 
ready used for standard CC could be used, but we will not investigate this further in the 
present article. 

Suppose the ground state is desired as initial condition. One option is to perform a 
Hartree-Fock calculation to generate a set of orbitals, and then perform a CCD calculation 
within this basis. This is the standard practice for molecular calculations, and even though 
not an exact critical point of the OATDCC energy, it should be a suitable starting point 
for dynamics calculations. In our numerical experiment in Section I VI I we choose a similar 
approach. 

VI. A NUMERICAL EXPERIMENT 
A. Outline and model system 

A numerical experiment on a model system mimicking electron-atom collision has been 
performed in order to test the OATDCCD method against the standard MCTDHF method. 
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(Recall that OATDCCD is an approximation to MCTDHF.) Generic MCTDHF and OATD- 
CCD codes has been written from scratch. The MCTDHF code was tested against numerical 



experiments reported in the literature 3| and found to agree perfectly with these. The OAT- 
DCCD code uses the algebraic expressions listed in|A]for the equations of motion. The code 
is tested against MCTDHF computations for special combinations of L and N where the 
two ansatze are equivalent. Perfect agreement was found, indicating the correctness of the 
implement at ion . 

The numerical experiment consists of two phases: (1) preparation of the initial wave- 
function, and (2) propagation of this state from t — to t — tfi na i using both MCTDHF 
and OATDCC while monitoring some observables. In particular the energy Eh should be 
conserved at all times. 

The test system is defined as follows. Consider a model consisting of N electrons in one 
spatial dimension. The orbitals are functions ip p (x,s), where x G K is the spatial position 
and s E { — is the quantum number of the projection of the electron spin along some 
arbitrary axis. 

The particles interact via a smoothed Coulomb force (for simplicity) and an external 
Gaussian well potential. The Hamiltonian of the system has one-body part 

h = -~ + V(x) t V(x) = -V e-^\ 

where Vq — 7 and a = 1.5 are the parameters of the Gaussian well. The smoothed Coulomb 
interaction is given by 

u(xx,x 2 ) 



V|S1-S2|+* 2 " 

and we use parameters A = 1 and 5 = 0.2. These are reasonable parameters for, say, a 
quantum wire model 



The complete Hamiltonian is seen to commute with any spin operator. We remark that 
if (1) the initial orbitals are spin-orbitals on the form ip p (x,s) = ip p (x)xcr p (s) where \ G is 
a spinor basis function, and if (2) the initial wavefunction is an eigenfunction for the total 
spin projection operator, the equations of motion ( 13 Op preserve these properties. That is to 
say, under these conditions the orbitals are always on product form 

tp p (x,s,t) = ip p (x,s,t)xa p (s), 

and the wavefunction remains an eigenfunction of the total spin projection. 
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B. Discretization and propagation scheme 



We discretize the one-particle coordinates by introducing a standard discrete Fourier 
transform-based discretization over the interval [—R,R] using iVgHd = 64 points {37]. The 
total number of basis functions is then Nb = 2N gT id = 128. The kinetic energy operator is 
evaluated using the fast Fourier transform (FFT) and is highly efficient and accurate. In 
our calculations, we set R = 15. The orbital matrices $ and <3> become standard matrices of 
dimension L x Nf, and 7V& x L, respectively, giving a simple representation in the computer 
code. 

For propagation, we choose a variational splitting scheme 



Variational splitting is a 



generalization of the standard split-step scheme used for brute-force grid discretizations of 
few-body problems [37 1. Using this scheme, the time step At can be chosen relatively large 
and independently of the grid spacing Ax. Variational splitting is most easily described in 
terms of a time- dependent Hamiltonian H(t) = T + 'Y^ = _ 00 5(t — (n + \/2)At)(H — T), where 
T is the kinetic energy operator and H— T is the remaining potential terms. A time step then 
consists of three steps: (1) Propagation using kinetic energy only as a time step At/2, (2) 
Propagation using H—T only a time step At using a fourth order Runge-Kutta for simplicity, 
and finally step (1) is repeated. The key point is that the stability of the scheme becomes 
insensitive to Ax, and that step (1) is evaluated exactly without involving the amplitudes 
at all. The local error is of order At 3 . A simple integration of Eqns. (|30|) using, say, Runge- 
Kutta will require a time step At ~ Ax 2 . In multi-configurational time- dependent Hartree 
calculations, other ways of eliminating this stability problem than variational splitting is 
often used, e.g., the constant mean-field scheme 2j. However, this requires more coding 
effort than variational splitting which is sufficient for our modest purposes. 



C. Preparation of initial wavefunction 

We need initial wavefunctions for both the MCTDHF and OATDCCD ansatze. The 
former is computed as follows. 

The parameters of the Gaussian well potential and interaction potential are experimen- 
tally chosen so that they support an N = 4 ground state ^4) with total spin projection zero. 
This ground state is computed in the MCTDHF scheme using imaginary time propagation 
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FIG. 2. Particle density n{x) of initial wavefunction and Gaussian well V(x). Units are arbitrary. 

of the equations of motion from a random initial condition. 

In order to generate non-trivial but easily understood dynamics, we prepare a fifth particle 
in a classical-like Gaussian wavepacket g(x, s) given by 

g(x, s) =C exp [-(x - x ) 2 / (4a 2 ) + ik Q x] X+1/2O), 

where C is a normalization constant, x the "starting position" of the particle, and ko is 
the "starting momentum". The parameter a controls the width of the wavepacket. For our 
experiment, we choose xq = 10, ko = 1.2, and, a = 1.25. Note that we (arbitrarily) choose 
spin +1/2 for this particle. The initial MCTDHF state is then 

|*>=^|*4), 

where is the creation operator associated with g(x, s). The complete wavefunction there- 
fore intuitively describes an incoming electron on collision course with a bound beryllium-like 
"atom" in the ground state. 

In terms of the MCTDHF parameters, the addition of the particle simply corresponds to 
extending the orbital matrix $ with an extra column (orthogonalized against the others), 
in addition to a lot of zeroes in the MCTDHF coefficient vector A. 

The spatial particle density n(x) = n(x, —l/2)+n(x, +1/2) of |\&) is shown in Fig.[2]along 
with the confining potential. The density n(x, s) is given by the diagonal of the reduced 
one-body density matrix 7(2:, s, x', s'), viz, 

j(x,s,x',s') = (^\^(x, s)il>(x',s')\^}, 
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or in terms of the coefficients p^ = [p^] and the orbitals, 

7 = (32) 

Having obtained the MCTDHF initial wavefunction, the OATDCCD initial condition is 
computed as follows. We first transform the orthonormal orbitals $ to generate orthonormal 
so-called Brueckner orbitals $b- The wavefunction amplitudes are transformed accordingly. 



By definition [39(j , the Brueckner orbitals optimize the overlap with the reference determi- 
nant, which is actually equivalent to the convenient property that the singles amplitudes 
vanish identically, the so-called Brillouin-Brueckner theorem. We now have 

|tf> = (0 b |W)(1 + A 2 + A 3 + ---)|0b) = (0 B |^>e T2+T3+ -|<M, (33) 

where |0b) is the determinant that has maximum overlap with \^f). We now simply take the 
OATDCC orbitals to be $b (and $ = $b), and let T 2 be as in Eqn. (I33I) . i.e., we perform 
a projection in CC amplitude space. 

It remains to define A 2 . We note that for the MCTDHF limit of OATDCC, = 
(^|/(^|^), which gives A n such that 

(^|e T ^ +r 3+- = (0 B |(1 + Ai + A 2 + ■••), 

which is used to extract A 2 algebraically in the computer code. 

We comment that the initial condition computed in this way only is a critical point of the 
coupled cluster energy to within an approximation, albeit a very good one. Alternatively, 
we could solve the CCD equations in the Brueckner basis, which gives very similar results. 

We comment that the OATDCC counterpart of the density matrix ( 1321) is 



7 



(34) 



The density plot of the OATDCC initial condition is visually indistinguishable from the 
MCTDHF initial condition, so we do not plot it separately in Fig. [2j 

D. Results 

Having obtained initial conditions, these are propagated in time with At = 0.005 until 
t = iflnai = 30. The energy is conserved to a very high precision, see Fig. [3j The energy- 
conservation also improves with reduced time step, as is expected. 
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FIG. 3. (Left) Energy conservation in the MCTDHF and OATDCCD schemes for the numerical 
experiment. The deviation from the initial energy is shown. The initial energies are -Emctdhf(O) = 
— 12.2102145 and -Eoatdccd(O) = —12.2104173. The gap between these numbers are due to 
neglection of the triples amplitudes and higher in the ground state. (Right) Integral of imaginary 
part of particle density, see Eqn. (j3"5|) . 

As discussed in Section III Bl expectation values may gain small imaginary parts. Along 
the CC computation we therefore monitor 

f( t ) = Y, / \lmn(x,s,t)\dx. (35) 

s 

In Fig. E3/(i) is displayed, and it is indeed a small number compared to the particle number 
N = 5. 

In Fig. H] we show the density as function of t of each calculation side by side. They are 
seen to agree qualitatively. The density evolution clearly shows how the incident electron 
interacts with the beryllium "atom" . Some of the density is clearly transmitted and reflected 
from the atom, while the atom is slightly perturbed, performing small amplitude oscillations. 
This demonstrates that manybody effects in the simulation are significant, and that the 
OATDCCD calculation captures these well. 

Quantitatively, the densities show some differences after the collision event that may look 
like a phase shift in the oscillation of the atom part. We have not investigated further, but 
conjecture this to be a result of the fact that the CC approximation changes the spectrum 
slightly. The absolute value of the density difference is shown in Fig. |5j 
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FIG. 4. Electron density plot as function of time for each calculation. MCTDHF is on the left, 
while OATDCC is on the right. The x-axis is divided into spin up (left half) and spin down (right 
half). The two densities are seen to be very similar, see also Fig. [5j The incident electron is 
clearly reflected and partially transmitted through the initially stationary beryllium atom. After 
the collision, the atom is seen to exhibit oscillations. The interference fringes at the end of the 
simulations are due to boundary effects. 

VII. CONCLUSION 

The bivariational principle for Schrodinger dynamics has been discussed at length, and the 
orbital-adaptive time-dependent coupled-cluster method (OATDCC) was developed. The 
method can be viewed as a systematic hierarchy of approximations to the highly successful 
multiconfigurational time- dependent Hartree method for fermions (MCTDHF), with simple 
time-dependent Hartree-Fock as the simplest case. The doubles approximation (OATD- 
CCD) was discussed in detail, and numerical experiment performed showing that the method 
gives sensible results. OATDCC scales polynomially where MCTDHF scales exponentially 
with the number of particles N. 

It was observed that imaginary time propagation for eigenvalue computation does not 
seem feasible for the OATDCC method. Studying methods for solving the time-independent 
orbital- adaptive CC should be a useful line of research. Such an eigenvalue computation 
method would constitute a hierarchy of approximations to the multiconfigurational Hartree- 
Fock method. 

OATDCC is easily generalized to bosonic systems, with particularly interesting appli- 
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FIG. 5. Differences between the OATDCC and MCTDHF electron densities calculated. The main 
feature is a phase-shift in the atom oscillations. 

cations to Bose-Einstein condensates (BEC). The resulting method will approximate the 



MCTDH method for bosons 40], and can possibly treat substantially more particles. In 
fact, the CC approximation is very well suited to describe a BEC since it naturally captures 
the idea of excitations on top of a condensate. 
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Appendix A: Algebraic expressions for CCD 

Here we list algebraic expressions for various quantities appearing in the OATDCC 
method using a doubles only ansatz, i.e., CCD. The expressions are computed using the 



second quantization toolbox in the Python library SymPy 
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22]. 



The expression for the expectation functional becomes 

_d_ 

\He T \<t>) + U4>»\e~ T He T \<P) (A2) 



S H [X, r, $, $] = (0|#e T |0> + ]T A M -^^[A, r, $, $] (Al) 



JHf\4>) + \^%{i$\e- T He T \<l>) (A3) 

where we have used linearity of £jj in A M , and where the latter expression explicitly states 
the expansion in the CCD case. 

To obtain computational formulae, we consider separately the one and two-body terms 
in H, i.e., 

H^ = J2h p q 4^ K = {Cp p \HV\y q ) (A4) 

pq 

and 

Hi2) -\j2<A c ic s c q , u v q : = (^ r \H {2) \^ s )As (as) 

pqrs 

The coefficients u^ r s are the the anti-symmetrized two-body integrals, see Equations ( TT2i) and 
(H3J). We get 

(0|# {1) e T |0> = /i{ (A6a) 
(0l^ {2) e T |0) = ^^ + l4 (A6b) 

for the CCD energy, and 

= ~h a A c P(ab) + fc 6 P(i?) (A7a) 



+ T%u bk P{ab) - r^tu k iP{ab) + r%u? c P(ab)P(ij) + \r%ut + < (A7b) 



for the derivatives with respect to X^ b . Inserting these expressions back into Eqn. flAip . the 
complete expression for the CCD expectation value functional becomes 

£ H [X, r, <£, $] = ^A£t£ + \H>%>t% - \xLt£u% + ^T^T^ig (A8) 



16 6 J ' c 8 iJ- 8 kj dc + 2 ab v ck 



4 

+ ^51<«£ + ^A>« 6 + ^fu* + £«g (All) 
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To solve the equations of motion for the amplitudes A^ 6 we also need the derivatives of £h 



ab 

. ^ Pnr thp nnp-VinHv nart 





with respect to r°f '. For the one-body part 



ab £ HW = ti k \itP(ij) - K\gP(ab), (A12a) 

ij 



and for the two-body part 
-£ 



ij 



+ l>&£ + l^ulP(ij) + \%u k \P(ij) - K h c r^ d P{ab)P{ij) 



+ ^ c < k P(ab)P(ij) + \\tri^ ah P{i 3 ) + \\^u% + iAfrig 
~ l&M d P(ab) + ul (A12b) 

The operator P(ij) is an anti-symmetrizer: f(ij)P(ij) = f(ij) — /0 Z )> an d similarly for 
P(ab). The appearance of a P(ij) or a P(ab) should be ignored for the invocation of the 
summation convention. 

Our calculations for expectation values and derivatives are of course valid for any one- 
or two-body operator. 

Expressions for the reduced density matrices p q v are readily obtained from £ H m by using 

K = 6 p' 6 1>- We § et Pi = Pa = °> While 

A = ^ = 6*5l - ~5<«r fc t = 51- -A^t (A13a) 

Pa = £ cUb = \&*>&$ = \^j\j (A13b) 

Finally, we compute the two-body reduced density matrix. We only list nonzero elements. 

p% = P(ij)*f*J - P(ij)P(kl)h^T^ + \\ k ^ (Al4a) 



2 t ca j m ' 2 



+ P(ij)l^rf k + \x k c l Air-f + 4 ( Al4b ) 
Pt = -Pt = -A = P«\ = \H>&% ~ K k c r% (Al4c) 



Pi = ^ (A14d) 
1 

2' 



PS* = o W (A14e) 
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Appendix B: Derivation of MCTDHF 



The purpose of this section is to provide a brief derivation of the MCTDHF method inde- 
pendently of the existing derivations in the literature, using the time-dependent bivariational 
principle. 

The bivariational MCTDHF manifold is 

M MCTDHF — |J V[<£] x V[$], (Bl) 

where the union is taken over all biorthogonal choices of orbitals. That is to say, is an 
arbitrary vector in the discrete (FCI) Hilbert space generated by $, and is an arbitrary 
vector in the discrete space generated by the dual orbitals $: 

|*> = 5>^> (B2) 
(§| = ]Ti,(<A (B3) 

where we include the reference determinants in the sums by defining |0 O ) = 10) and (0o| = 
i.e., fi = is the reference. Since the FCI spaces are linear it is allowed to remove the 
phase and normalization ambiguity in Eqn. (j3J) by requiring (^l^) = 1. 
The time derivative of is 

J^) = ^>|0 M > + iW, (B4) 

where D is defined in Eqn. (TPTI) . We obtain two expressions for the action functional, which 
are equivalent: 

S[A, A, $, $] = f i/L4 M > + \%A il D^A v - A fl H^A u dt (B5) 
Jo 

\KA^ + pl(ihrf q - hP q ) - \pi;y q : dt, (B6) 



where 



and where 



D^ = {r\D\cf> v ), H^^HW^), r£ = (0 p \ip q ) t (B7) 



Pl = (*|^|*>, p q p s r = (*|ct c tc s c 9 |*> (B8) 
are the reduced one- and two-body density matrices. 
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In a similar way as for the CC ansatz, it can be shown that mixing of the orbitals $ — > $G 
and $ — > G~ x &, G being an arbitrary invertible L x L matrix, with a corresponding inverse 
operation on A and A leaves the action invariant (except for a total time derivative with 
vanishing variation). This implies that all the inner products if can be chosen arbitrarily, 
and we choose ri v q = for simplicity, which implies D^; = 0. This corresponds to the usual 
gauge choice in MCTDHF theory, see Ref. The resulting permissible variations in (p p 
are then of the form 5ip p = x — QXi with Q — 1 — $$. Similarly, the permissible variations 
in (p p are 5<p p = x = xQ- The P-space variations are identically zero in this gauge. 

Performing first the variations in A v and A M , respectively, we obtain the equations 

\KA V = H»A 11 and -ihX^ = A u H», (B9) 

i.e., the time-dependent Schrodinger equation in a moving basis and its dual equation. Since 
H v = ((j) l '\TlHTl\(j) l j), the operator expression (llip is useful for implementations. 

Performing the variations in (p p we obtain, in a manner similar to Section IIVBI the 
equation 

q 

and performing the variation in (p q we obtain 

p 

which are identical in form to Eqns. (I30ep and (I30fl) . 

Assume now that $(0) = $(0) H at time t = 0, and that A M (0) = A^(0)*. Then clearly 
Pl = (p?)* and Pp % = (ff q :y. Similarly, h\ = and u% = (u£)*. 

From this it follows that (<p p \ = \<p p )\ \ = (A' 1 )* and hence (&(t)\ = |\E r (£)) t for all t. 

Comparing with the equations of motion in, say, Ref. [4^], we see that we indeed have 
arrived at the MCTDHF equations of motion. 

Appendix C: Some techinical proofs 

1. Orbital transformations 

We consider the orbital transformation 

$ — ► $' = $G, $ — ► = G- 1 ®, (CI) 



q qrs 



;bio) 



^2p q p((p P \h+^2p pr 



prs 



Q. 



(Bll) 
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where G = exp(g) is an L x L invertible matrix. (Recall that any invertible matrix has a 
logarithm.) The transformation preserves biorthogonality. The transformation is equivalent 
to transforming the creation and annihilation operators as 

We wish to show, that Eqn. ( 1C2|) is again equivalent to 

4 — >• e§ 4 e ~^ c g — ► e% e - § , (C3) 

where 

# = E ^ P cJ £ p- ( C4 ) 

We begin by noting that the similarity transform can be expanded using the BCH formula 

as 



1 «— ■ 1 

e 9 xe~ 9 = x+[g,x} + —[g,[g,x]] + --- = ^2— [g, x) n (C5) 

n=0 

where [A, B] n is the n-fold nested commutator. Suppose we can show, that for all n > 0, 



and 



&,p I ]» = X)(- 1 ww ( C6b ) 



This would imply 



e 9 c\e~ 9 

" n 



oo 1 

E^E(A P ^E^ c ! ( C7 ) 



and 



n=0 <j 



e%e~ 9 = K —T 5> n )<^ = E ( C * 



ft! 

n=0 p p 



which would prove the result. 

We prove Eqn. ( 1C6|) by induction, and for simplicity we consider only ( lC6aj) . The proof 
for (]C6b|) is similar. For n — 0, = 5 gp , which shows that n = is trivially true. 
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Suppose now, that Eqn. (!C6aP holds for some n > 0. We show that it holds for n + 1 as 
well. We get 

[&cj] n+ i = [g, [g,c\] n ] = lg,^2gr q ct] = ^g sp gr q [c%,4\ = ^^^4 = ^2dsp lc l 



psr 



psr 



which proves the induction step. We have used that 

which is easily calculated using the fundamental anticommutator. 



(C9) 



(CIO) 



2. Amplitude equations 



In this section, we consider the derivation of the amplitude equations ( EM|) from the 
variation of the OATDCC functional with respect to the amplitudes A M and r M , which are 
all independent variables. We therefore start with varying a single amplitude A„(t), i.e., 
SX^t) = for /i^y, and 5r^{t) = for all \i £ X. Note that all variations vanish at the 
end points t = and t — T of the time interval. We wish to compute 



SS[\, t}=5 [ ih\„{t)t»{t) - Sh-vkdoW), r{t)\ 
Jo 



(Cll) 



where we have suppressed the dependence on $ and $ in £ and S since they are held fixed 
in the variation. We obtain 



SS[X, r] 



T ih5X v (t)f v - dW^(t),r(t)] ^ (t) ^ 



.... d£ H [X(t),r(t)} 



dt = 0. 



(C12) 
(C13) 



Note that we do not invoke the summation convention. Since the expression must vanish 
for all choices of X u (t), it follows that 



m v {t) 



_d_ 

dX, 



-~ihD a 



(C14) 



must hold. 

Similarly, we hold A^ fixed, and vary only a single amplitude t 1 



SS[X, t] 



T 



5r v {t) 



-ihX,, 



d£ H -ihD [Ht),r(t)] 
dr v 



dt + ih[X u (t)Sr u (t)} 



T 

t=0 



(C15) 
(C16) 



47 



where we have performed integration by parts. The boundary terms vanish by definition of 
the variation, so that we get 

-mx v {t) = ^£ H -w%[\(t),T(t)Mt)Mt)]- (Ci7) 
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